
clear all
************************************
*Numeric Export disabled
local TEX ""
local TEXON 1 //Change "0" to "1" to enable exporting numeric results to .txt files
if !`TEXON' local TEX "*"
*Also need to install `texresults2
*ssc install texresults2
************************************

*************************
*set working directory
*************************

clear all 

use "Enns_SociologicalScience_Reproduction\DS1_0.dta", clear 

set matsize 10000

gen policyadopted =.
recode policyadopted .=0 if OUTCOME==0
recode policyadopted .=1 if OUTCOME==2
recode policyadopted .=1 if OUTCOME==3
recode policyadopted .=1 if OUTCOME==4


*****************
gen  pred50_sw_lor = log(pred50_sw/(1- pred50_sw))
label var pred50_sw_lor "log of odds: pred50_sw"

gen  pred90_sw_lor = log(pred90_sw/(1- pred90_sw))
label var pred90_sw_lor "log of odds: pred90_sw"

gen  pred10_sw_lor = log(pred10_sw/(1- pred10_sw))
label var pred10_sw_lor "log of odds: pred10_sw"
*******************


***********************************************************
//Reproduce Gilens' multivariable results (NOT DEFLATED) Table A3.3 P.256, to confirm approach used in simulations.
drop if policyadopted==.
*Save DV as vector
mkmat policyadopted, matrix(y)

drop if pred10_sw==.
drop if pred50_sw==. 
drop if pred90_sw==.

gen cons=1
mkmat cons pred10_sw pred50_sw pred90_sw, matrix(X)

*Need to standardize so mean = 0
foreach var of varlist pred10_sw pred50_sw pred90_sw {
sum `var'
gen z`var' = `var'-`r(mean)'
}

mkmat cons zpred10_sw zpred50_sw zpred90_sw, matrix(Xz)

mat def XpXz = Xz'*Xz

mat def b = invsym(XpXz)*Xz'*y
*TABLE A3.3 P.256 (Column 2)
mat l b
***********************************************************


***********************************************************
*Repoduce deflated results, Table A3.3 P.256 (Column 3) to confirm simulation approach
mat deflate = (1, 1, 1, 1 \ 1, 1, .86, .81 \ 1, .86, 1, .83 \ 1, .81, .83, 1)
*gen empty matrix
matrix XpXzdef = J(4,4,0)

forvalues i = 1/4 {
	forvalues j = 1/4 {
		matrix XpXzdef[`i',`j'] = XpXz[`i',`j']*deflate[`i',`j']
		}
	}
mat def b = invsym(XpXzdef)*Xz'*y
mat l b
***********************************************************


******************************************************************
******************************************************************
// Appendix 10.2: **How sensitive are Gilens' multivariabel resuls to the specific estimate of correlated measuremeent error. This is important, becuase Appendix 10.1 shows the estimate of correlated measurement error is necessarily inflated. ment error
******************************************************************
******************************************************************
//save estimates
tempfile sims_multivariate3
tempname temp1
postfile `temp1' def1090 def5090 def1050 ///
			 multi10 multi50 multi90 ///
			multi10se multi50se multi90se ///
				using "`sims_multivariate3'", replace
forval v = 1/18 {
scalar def1090 = .81+`v'/100
forval v2 = 1/16 {
scalar def5090 = .83+`v2'/100
forval v3 = 1/13 {
scalar def1050 = .86+`v3'/100
mat deflate = (1, 1, 1, 1 \ 1, 1, def1050, def1090 \ 1, def1050, 1, def5090 \ 1, def1090, def5090, 1)
*gen empty matrix
matrix XpXzdef = J(4,4,0)

forvalues i = 1/4 {
	forvalues j = 1/4 {
		matrix XpXzdef[`i',`j'] = XpXz[`i',`j']*deflate[`i',`j']
		}
	}

mat def b = invsym(XpXzdef)*Xz'*y
mat l b
local multi10 = b[2,1]
local multi50 = b[3,1]
local multi90 = b[4,1]
local multicons = b[1,1]

mat def u = y-Xz*b

mat def sumsq = u'*u
scalar sumsq = sumsq[1,1]
local sumsq = sumsq
di `sumsq'
local sigma2 = `sumsq'/(1779-4)
di `sigma2'

mat def varcov = `sigma2'*invsym(XpXzdef)

matrix se = J(4,1,0)

forvalues i = 1/4 {
		matrix se[`i',1] = sqrt(varcov[`i',`i'])
		}
mat l se
local multi10se = se[2,1]
local multi50se = se[3,1]
local multi90se = se[4,1]
local multiconsse = se[1,1]
//post results to dataset 
post `temp1' (def1090) (def5090) (def1050) ///
			(`multi10') (`multi50')  (`multi90') ///
			(`multi10se') (`multi50se') (`multi90se')
}
}
}

postclose `temp1' //write results to dataset 

********************************************
use "`sims_multivariate3.dta'", clear
********************************************

gen t50 = multi50/multi50se
gen t10 = multi10/multi10se
gen t90 = multi90/multi90se

count if t50>1.96 | t10>1.96
local pos = r(N)
count if t10<-1.96 | t50<-1.96 | t10>1.96 | t50>1.96
local imp = r(N)
count
local tot = r(N)

local percpositive = 100*`pos'/`tot'
local percimplausible =  100*`imp'/`tot'

`TEX'texresults2 using TxtFiles_NumericalResults\multivariablesims.txt, texmacro(percpositive) result(`percpositive') round(1) replace //
`TEX'texresults2 using TxtFiles_NumericalResults\multivariablesims.txt, texmacro(percimplausible) result(`percimplausible') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\multivariablesims.txt, texmacro(totaln) result(`tot') round(0) append
